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Abstract 

The  isothermal  fuel  cell  model,  given  in  an  earlier  publication,  will  be  generalized  to  describe  the  behaviour  of  non-isothermal  fuel 
cells  of  co-flow  type.  To  this  end  the  temperalure  distribution  inside  a  fuel  cell  in  steady  state  is  investigated  analytically.  A  simplified 
relation  between  the  local  temperature  and  the  fuel  utilization  is  derived  and  its  practical  significance  elucidated.  Furthermore,  it  is  shown 
that  the  solution  of  the  non-isothermal  model  is  accurately  approximated  by  analytical  expressions  obtained  from  a  so-called 
quasi-isothermal  approach.  This  new  approach  yields  a  similar  expression  for  the  cell  voltage  as  derived  from  the  isothermal  model.  The 
quasi-isothermal  approach  is  also  used  to  make  a  clear  comparison  between  the  isothermal  and  the  non-isothermal  fuel  cell  model.  ©  1998 
Elsevier  Science  S.A. 
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1.  Introduction 

In  a  previous  article  [1]  the  setup  of  an  isothermal  model  for  a  fuel  cell  in  steady  state  was  given.  In  this  model,  the  local 
fuel  utilization  u  and  the  cell  voltage  Vtcll  are  the  solution  of  an  appropriate  ‘boundary  value  problem’  (BVP).  based  on  the 
following  quasi-ohmic  relation  for  the  local  current  density  i 

Kq(«)  “  Fell  =  '''(«)  (0 

where  the  local  equilibrium  potential  V,,Ju)  depends  on  position  through  the  fuel  utilization,  u,  whereas  the  quasi-ohmic 
resistance,  r,  was  assumed  to  be  constant. 

An  analytical  approach  for  isothermal  fuel  cells  was  introduced,  based  on  the  linearized  Nernst  equation 

K;(o)  =  Ve;(0)  -au  (2) 

The  relation  between  the  parameter  a  in  Eq.  (2)  and  the  Nernst  loss  was  shown  and  the  accuracy  of  the  analytical  results 
was  verified  by  comparison  with  numerical  computations  and  measurements  on  a  bench-scale  molten  carbonate  fuel  cell 
(MCFC;  1000  cm2). 

The  isothermal  model  can  be  applied  if  the  influence  of  temperature  gradients  on  the  local  equilibrium  potential  and  the 
quasi-ohmic  resistance  is  negligible.  The  local  equilibrium  potential  depends  only  slightly  on  temperature  [2,3],  but  the 
condition  of  a  uniform  quasi-ohmic  resistance  is  a  major  restriction  of  the  isothermal  model. 

The  objective  of  this  article  is: 

1.  to  develop  a  non-isothermal,  analytical  model  for  fuel  cells  with  co-flow,  and 

2.  to  compare  the  non-isothermal  model  with  the  isothermal  model. 

The  first  item  will  be  realized  by  generalization  of  the  isothermal  cell  model  to  a  quasi-isothermal  model.  The  setup  of 
the  quasi-isothermal  model  is  based  on  a  mathematical  method  that  is  similar  to  ‘Duhamcl’s  principle’  [4]  or  the  method  of 
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"variation  of  parameters’  [5].  The  quasi-isothermal  model  will  provide  similar  expressions  for  the  cell  voltage  of  a 
non-isothermal  fuel  cell,  as  obtained  ii  Ref.  [1]  for  an  isothermal  fuel  cell. 


2.  Non-isothermal  model  for  fuel  cellls  with  co-flow 

2.1.  Problem  definition 

As  in  the  isothermal  model  [1],  a  rectangular  flat  plate  fuel  cell  is  considered,  but  the  cell  temperature  and  the 
quasi-ohmic  resistance,  r,  are  no  longer  assumed  to  be  uniform.  The  other  assumptions  that  underlie  the  isothermal  model 
will  be  maintained  and  are  listed  in  Table  1.  The  cell  part  between  two  cross  sections  of  the  cell  is  called  ‘subcell’.  An 
example  of  a  subcell  is  shown  in  Fig.  I.  Since  the  current  collectors  will  not  affect  the  problem,  they  are  omitted  in  this 
figure.  Furthermore,  it  is  assumed  that  the  whole  cell  and  all  subcells  are  thermally  insulated.  During  the  setup  of  the  model, 
the  generality  of  this  latter  assumption  with  respect  to  insulated  fuel  cell  stacks  will  be  made  clear. 

For  a  fuel  cell  operating  in  steady  state  with  co-flow,  the  temperature  will  be  almost  uniform  over  a  cross  section  of  the 
cell.  For  instance,  temperature  differences  over  a  cross  section  of  an  MCFC  or  solid  oxide  fuel  cell  (SOFC),  with  co-flow, 
are  in  the  order  of  10  K  [6],  It  is,  therefore,  convenient  to  assume  a  uniform  gas  temperature  over  a  cross  section  of  the  cell 
(assumption  7  in  Table  2).  Denoting  the  local  temperatures  of  the  anode  and  cathode  gases  by  Ta  and  Tc,  respectively,  this 
means  that 

T(x)=Ux)=Te(x) 

Since  Tfix)  and  Tc(x)  are  not  really  identical,  the  proper  Tlx)  actually  represents  a  mean  gas  temperature,  Nevertheless, 
the  local  equilibrium  potential  depends  only  slightly  on  temperature  and  will  be  a  function  of  T(x)  and,  of  course,  a 
function  of  the  local  gas  compositions,  in  a  very  good  approximation.  The  quasi-ohmic  resistance  is  usually  very  sensitive 
to  temperature,  but  in  a  first  approach  it  can  also  be  considered  as  a  function  of  T.  Moreover,  the  expression  that  will  be 
derived  for  the  local  gas  temperature,  T,  requires  only  a  slight  modification  in  order  to  approximate  the  hardware 
temperature.  Hence,  the  small  loss  in  accuracy  by  considering  the  quasi-ohmic  resistance  as  a  function  of  T  is  easy  to 
compensate. 

The  local  fuel  utilization  in  the  non-isothermal  model  will  be  distinguished  from  the  corresponding  variable  in  the 
isothermal  model,  by  w'riting  u  and  u,  respectively.  As  in  the  isothermal  model,  the  local  gas  compositions  can  be 
calculated  from  the  local  fuel  utilizat  on  it.  So,  neglecting  the  hydrostatic  pressure  gradients  (assumption  6  in  Table  1),  the 
local  equilibrium  potential  Veq  and  the  local  quasi-ohmic  resistance  r  can  be  considered  as  functions  of  u  and  7. 
Analogously  to  Eq.  (1),  the  expression  for  the  local  overpotential  can  still  be  written  in  accordance  with  Ohm  s  law',  as 

Kq[SJ]-Vcell  =  it«,r]r[n.T]  (3) 

With  the  use  of  this  current-voltage  relation,  the  BVP  derived  in  Ref.  [1]  can  be  generalized  to  describe  also  the 
non-isothermal  behaviour  of  fuel  cells. 

It  is  convenient  to  repeat  the  original  form  of  this  BVP  first 

~  =  —{Kq[u  Jal[]~Vccl]  ),  (0  <  x  <  1) 

dx  n,„ 

!<(0)  =0,  u(  1)  =  U  j 

with  x  the  scaled  distance  to  the  cell  inlet,  it,  the  total  fuel  utilization,  iin  the  equivalent  input  current  density  of  the  fuel 
supplied  and  Tcen  the  uniform  cell  temperature. 


BVP(  1)  (isothermal) 


Table  1 

Overview  of  the  basic  mode!  assumptions,  equal  to  those  for  the  isothermal  model  [1] 

1  Stationary  fuel  cell  (time  independent) 

2  Large  oxidant  flow;  negligible  oxidant  utilization  with  respect  to  the  Nernst  loss 

3  Changes  in  fuel  composition  are  significant  in  the  direction  of  the  cell  outlet  only 

4  The  local  overpotential  depends  linearly  on  the  local  current  density  through  a  quasi-ohmic  resistance 

5  All  reactions  in  the  gas  phase  are  in  equilibrium 

6  Negligible  pressure  drops  (uniform  pressure) 
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Fig.  1.  Longitudinal  section  of  the  cell.  Without  loss  of  generality  the  setup  of  the  cell  model  is  restricted  to  an  insulated  cell  for  which  the  heat  transfer  Q 
is  equal  to  zero. 


Analogously  to  the  derivation  of  BVPfl),  it  can  be  shown  from  Eq.  (3)  that  the  local  fuel  utilization  u{ x)  and  the  cell 
voltage  Vcell  now  satisfy 


d  u 
d  a 


5(0)  =0,  h(  1)  =  u  j 


(0<jc<  1) 

/BVP(2)  (non  —  isothermal) 


In  BVP(2)  the  local  temperature  and  the  quasi-ohmic  resistance  depend  on  position,  whereas  in  BVP(1)  these  quantities 
are  uniform.  Evidently,  for  arbitrary  quasi-ohmic  resistances,  r  and  r,  the  cell  voltages  in  BVP(l)  and  BVP(2)  are  not 
equal.  Therefore,  the  cell  voltage  in  BVP(2)  is  also  marked  with  a  tilde. 

When  solving  BVP(2),  we  will  consider  the  local  temperature  as  a  function  of  the  local  fuel  utilization  u,  rather  than  of 
the  position  x.  It  will  be  shown  that  the  description  of  T  in  terms  of  u  is  relatively  simple  and  very  effective.  Once  BVP(2) 
is  successively  elaborated  with  an  expression  for  T{u)  and  solved  for  the  set  (w(-v),Vcell),  the  local  current  density  i(x)  can 
be  calculated  from  the  same  identity  as  derived  for  the  isothermal  model 


d  it 

'( x)  =  (0  <-*<!) 


(4) 


In  a  fuel  cell  stack  the  operating  conditions  for  the  individual  cells  will  vary.  Nevertheless,  test  results  suggest  that  the 
influence  of  stacking  on  the  cell  performance  is  limited  by  applying  sophisticated  stack  designs  [7-9],  In  that  case  only  the 
performance  of  the  cells  at  top  and  bottom  of  the  stack  will  be  essentially  different  from  the  other  cells.  Consequently,  for 
an  insulated  stack  with  co-flow  it  suffices  to  consider  a  single  cell,  imposing  adiabatic  boundary  conditions  [6,10,1 1],  In 
fact,  if  the  net  heat  release  through  the  separator  plates  of  an  individual  cell  would  not  be  zero,  heat  would  be  released  from 
all  cells  (because  of  their  equivalence)  and.  hence,  from  the  insulated  stack  as  a  whole,  which  is  a  clear  contradiction. 
Consequently,  it  is  irrelevant  whether  the  separator  plates  inside  the  stack  are  thermally  insulating  or  conducting,  since  w'e 
use  a  uniform  (i.e.  mean)  temperature  over  a  cross  section  of  the  cell.  Though  the  setup  of  the  model  is  given  for  an 
insulated  cell,  it  also  holds  approximately  for  one  cell  of  an  insulated  stack.  Since  thermal  insulation  of  a  whole  stack  is 
subject  to  practical  difficulties,  many  literature  data  are  dedicated  to  numerical  stack  models  that  also  simulate  the  influence 
of  heat  transfer  from  the  stack  surface  to  the  surroundings  ([12-14]).  Furthermore,  in  Ref.  [15]  the  temperature  distribution 
in  a  fuel  cell  stack  of  co-flow  type  is  analytically  investigated  under  the  assumptions  of  a  fixed  temperature  at  the  stack 


Table  2 

Extension  of  the  basic  model  assumptions  to  the  non-isothermal  model,  and  additional  assumptions  (10  and  11)  that  underlie  the  linear  relation  between 
temperature  and  fuel  utilization 

7  Anode  and  cathode  gases  almost  have  the  same  local  temperature,  i.e.  the  gas  temperature,  T, 
is  considered  to  be  uniform  over  a  cross  section  of  the  cell 

8  Thermally  insulated  cell 

9  The  cell  surface  is  much  larger  than  a  cross  section  of  the  cell  (practically  insulated  .subcells) 

10  Negligible  heat  production  by  non-electrochemical  reactions  in  the  gas  phase 

1 1  The  total  heat  capacity  of  the  process  flows  is  almost  uniform  (large  oxidant  flow) 
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boundary  and  a  uniform  heat  generation  inside  the  stack.  Nevertheless,  modeling  results  show  that  the  major  impact  of  heat 
release  from  a  stack  is  restricted  to  the  cells  near  the  stack  boundary. 

In  general,  the  cell  surface  is  much  larger  than  a  cross  section  of  the  cell  (assumption  9  in  Table  2).  Because  of  the 
relatively  small  cross  section,  conduction  of  heat  will  hardly  contribute  to  the  total  transport  of  heat  in  the  flow  direction  of 
the  gases  '.  It  is  realistic  to  assume  that  transport  of  heat  in  the  flow  direction  happens  by  means  of  convection  only.  Then, 
all  subcells  may  be  considered  as  being  insulated  too.  It  is  noted  that  heat  conduction  perpendicular  to  the  flow  direction  is 
significant,  since  it  has  a  major  contribution  to  the  transfer  of  heat  from  the  hardw-are,  i.e.  the  electrode-electrolyte 
assembly,  to  the  gas  flows. 

Since  all  subcells  are  essentially  insulated,  the  local  increase  in  temperature  of  the  process  flows  can  be  calculated  from 
their  heat  capacities  and  the  local  heal  production.  In  Section  2.2  the  details  of  this  step  will  be  given  and  a  relation  for  the 
local  temperature  in  terms  of  u  will  be  introduced.  In  Section  3,  a  mapping  u  — »  n,  that  transforms  the  non-isothermal 
problem  BVP(2)  to  the  isothermal  problem  BVP(l).  will  be  introduced.  The  transformed  problem  will  be  referred  to  as  the 
quasi-isothermal  problem  and  can  be  used  to  generalize  the  results  and  conclusions  that  are  given  in  Ref.  [l]  for  an 
isothermal  fuel  cell. 

2.2.  Determination  of  the  local  temperature 

2.2.1.  A  general  expression  for  the  local  temperature 

Since  all  subcells  are  considered  to  be  thermally  insulated,  energy  is  transferred  across  the  border  of  a  subcell  by  the 
electric  current  or  the  gas  flows,  but  not  by  a  heat  flux.  The  entire  heat  production  inside  a  subcell  will  contribute  to  the 
increase  in  the  temperature  of  the  gas  flows.  Since  the  cell  is  operating  in  steady  state,  the  temperature  of  the  cell  hardware 
will  not  change.  Only  the  gases  are  heated  and  consequently  only  the  heat  capacities  of  the  gas  flows  will  be  relevant  to 
determine  the  temperature  rise  in  the  cell.  The  local  heat  capacity  of  the  anode  flow  is  denoted  by  mfu)cp  a(«),  where 
m.fu)  and  cp  fu)  are,  respectively,  the:  local  anode  mass  flow  and  the  corresponding  local  specific  heat  at  constant  pressure. 
Using  a  similar  notation  for  the  local  heat  capacity  of  the  cathode  flow,  the  total  heat  capacity  of  the  anode  and  cathode 
flows  is  given  by  mfu)cp  fu)  +  thfu)cp  fit).  This  heat  capacity  is  a  function  of  the  gas  compositions  only,  assuming  that 
the  anode  and  cathode  flows  are  mixtures  of  gases  with  constant  specific  heats  !.  It  determines  the  amount  of  heat  necessary 
to  achieve  an  increase  in  gas  temperature,  T,  of  one  Kelvin.  Hence,  the  amount  of  heat  related  to  the  temperature  rise  in  an 
insulated  subcell  at  the  cell  inlet  (see  Fig.  1),  is  given  by 

[mfu)cp/fu)  +mc(u)cp.c(u)][T(u)  -  T(0)]  ^  A //„[«, T- 7(0)]  (5) 

The  quantity  A//thermal  equals  the  amount  of  heat,  per  unit  of  time,  that  would  be  released  if  the  local  gas  flows  are 
cooled  from  the  local  temperature  7(5)  to  the  inlet  temperature  of  the  cell  7(0).  Below  we  will  derive  an  expression  for  this 
quantity,  from  the  first  law  of  thermodynamics. 

Neglecting  changes  in  potential  and  kinetic  energy  of  the  process  flow's,  the  first  law  states  that  the  heat  added  to  an  open 
system  in  steady  state  is  equal  to  the  change  in  enthalpy  AH  of  the  process  flows,  per  unit  time,  plus  the  power,  W.  that  is 
delivered  by  the  system  [17] 

Q  =  Houl-Hin+W^AH+W  (6) 

The  power  necessary  to  achieve  changes  in  the  volume  of  the  process  flows  is  taken  into  account  by  the  term  AH.  In  the 
case  of  a  fuel  cell,  the  term  W  represents  the  electrical  power  that  is  delivered  to  the  external  load  resistance,  i.e.  the  product 
of  the  cell  voltage  and  current.  Furthermore,  for  an  insulated  subcell  the  added  heat  per  unit  time  Q  is  zero.  So,  for  an 
adiabatic  subcell  at  the  cell  inlet,  in  which  a  fraction  u  of  the  fuel  supply  is  utilized,  the  first  law  is  expressed  by 

W(u)=  ~AH[it,T].  W(u)  =  VccUIinu  (7) 

Changes  in  thermodynamical  state  properties,  such  as  the  enthalpy,  are  independent  of  the  many  possible  ways  in  which 
the  change  of  state  can  be  achieved  [17],  For  example,  the  change  in  enthalpy  A H[7t,T]  could  equally  well  be  the  result  of 
an  isothermal  process  and  a  non-isothermal  process,  successively  causing  the  changes  in  gas  compositions  and  the  increase 
in  temperature.  In  Fig.  2,  these  subprccesses  are  schematized.  Without  loss  of  generality  we  may  state  that  the  work,  W,  is 


1  The  justification  of  this  statement  follows  from  a  dimension  analysis. 

2  The  specific  heats  at  constant  pressure  will  be  calculated  for  the  average  temperature  of  the  fuel  cell.  Changes  in  the  specific  heats  due  to  temperature 
rise  in  the  cell  will  be  neglected.  This  is  normally  allowed,  at  least  for  reasonably  limited  temperature  ranges  [16], 
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Fig.  2.  Schematic  representation  of  the  two  subprocesses  that  equal  the  change  in  enthalpy  of  the  gas  flows  through  an  insulated  subcell  at  the  cell  inlet. 


delivered  by  the  isothermal  subprocess  only.  Then  the  change  in  enthalpy  in  the  non -isothermal  subprocess  is  due  to  heating 
only.  In  other  words,  the  second  subprocess  yields  the  quantity  A  //[hcrmai  that  is  needed  in  Eq.  (5) 

A//[u,7-]  =  A//isothermal[«;r(0)]  +Aw,hermal[«,r-r(0)]  (8) 

In  Eq.  (8),  the  term  A//isothermal[«;7’(0)]  refers  to  the  isothermal  enthalpy  difference  due  to  changes  in  the  gas  compositions. 
The  quantity  -  A//isolhcrmal[«;7'(0)]  would  be  the  equivalent  of  the  released  heat  if  the  gases  would  be  converted  in  an 
isothermal  system  without  power  generation  (W  =  0).  However,  in  a  fuel  cell  the  change  in  enthalphy  is  controlled  by  the 
processes  in  the  electrode-electrolyte  assembly  and  -  A  Hhothfrmu,  is  partially  converted  into  electrical  power,  see  Fig.  2. 
Unfortunately,  it  is  not  fully  converted  into  electrical  power  and  the  remaining  part  still  represents  a  heat  production.  Since 
the  produced  heat  cannot  be  released  from  the  insulated  subcell  by  means  of  conduction,  the  temperature  of  the  gases  will 
increase.  This  gives  rise  to  an  extra  change  in  enthalpy,  which  is  represented  by  the  second  term  on  the  right-hand  side  of 
Eq.  (8).  The  first  law  (Eq.  (7))  relates  the  ret  change  in  enthalpy  to  the  delivered  power,  IV.  Substitution  of  Eq.  (8)  into  Eq. 
(7)  yields 


Vcell/In«  =  -  A/4 otherma,[«;7’(0)]  -  A //„[«, T-  7(0)]  (9) 

Using  Eq.  (9)  to  elaborate  Eq.  (5),  the  local  temperature  can  be  written  as 


T{ii:VccU)  =  T(  0)  + 


^isothermal  [«^(0)]  ~  ^cellAr." 


(10) 


Eq.  (10)  gives  the  local  temperature  inside  a  fuel  cell  as  a  function  of  the  local  fuel  utilization  u  and  the  cell  voltage,  if 
expressions  for  A//isotherma|[w,r(0)]  and  the  total  heat  capacity  of  the  process  flows  are  available.  In  the  next  section 
simplified  expressions  for  these  quantities  will  be  given. 


2.2.2.  A  simplified  expression  for  the  local  temperature 

The  obtained  expression  for  the  local  temperature  can  be  simplified  considerably,  if  heat  effects  due  to  chemical 
reactions  in  the  gas  phase  are  of  minor  importance  (assumption  10  in  Table  2).  For  instance,  in  general,  the  heat  effect  of 
the  water-gas  shift  reaction  is  small  in  comparison  with  the  reaction  heat  of  the  overall  cell  reaction.  Hence,  the  simplified 
expression  for  the  local  temperature  will  apply  to  fuel  cells  without  internal  reforming. 

Usually,  a  reaction  heat  at  constant  pressure  is  directly  proportional  to  the  amount  of  products  converted.  However,  the 
reaction  heat  A//lsothermal[tl;7(0)]  was  defined  with  respect  to  the  overall  process  in  the  gas  phase  and  is  not  related  to  one 
reaction.  Nevertheless,  in  the  case  that  the  heat  production  by  chemical  reactions  in  the  gas  phase  is  negligible  (assumption 


Table  3 


Standard  gas  compositions  in  mole  fractions  for  an  MCFC 


Anode  inlet  gas 

H  , 

CO 

h2o 

co2 

ch4 

Before  equilibrium 

0.64 

0 

0.20 

0.16 

0 

After  shift  (600°C) 

0.57 

0.07 

0.27 

0.09 

0 

Cathode  gas 

Air 

CO, 

Homogeneous 

0.70 

0.30 
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Table  4 

Parameter  values  for  an  MCFC.  corresponding  with  the  standard  gas  compositions  at  atmospheric  pressure 


X 

cpa(0) 

r,,(0) 

4, 

fna 

anode  flow  volume 

(V) 

(J/kg  K) 

(J/kg  K) 

(A) 

(kg/s) 

(cm3/s) 

1.299  a 

2935 b 

1 1 64 b 

1 

0.966X10  7 

0.59 

a  Value  corresponding  with  an  inlet  temperature  of  600  °C  and  a  total  fuel  utilization  of  95%.  Practical  values  of  y  depend  slightly  on  the  heat  effect  of  the 
shift  reaction.  Without  the  correction  made  for  the  shift  reaction.  y  would  be  equal  to  1.28  V.  irrespective  of  the  total  fuel  utilization  and  the  gas 
compositions  (see  Appendix  A). 

b  Values  determined  with  respect  to  the  average  cell  temperature  (650  °C).  The  values  for  y,  cp  ,  and  c  are  based  on  the  JANAF  tables. 


10  in  Table  2),  the  reaction  heat  ^  H,tol/iermal[u',T(0)\  will  essentially  refer  to  the  overall  cell  reaction  only.  Consequently,  it 
will  be  directly  proportional  to  the  fuel  utilization  u,  hence  to  the  corresponding  output  current 

A  isothermal  [“^(O)]  =  ~Xl  J'  (11) 

where  ^  is  a  constant  of  proportion  (with  dimension  V).  For  ideal  gases,  this  constant  depends  on  the  inlet  temperature,  and 
the  overall  cell  reaction  only.  In  Appendix  A,  the  value  of  x  >s  calculated  for  the  reaction  of  hydrogen  with  oxygen. 
However,  if  also  the  water-gas  shift  reaction  occurs,  assumption  10  is  not  exactly  satisfied  and  it  is  convenient  to  use  a 
slightly  modified  value  for  \  (see  Table  4  or  Appendix  A).  This  modified  value  depends  also  on  the  inlet  gas  compositions 
and  on  the  amount  of  converted  CO. 

Furthermore,  variations  in  the  total  heat  capacity  of  the  process  flows  will  be  small  in  the  case  of  a  large  and  almost 
homogeneous  oxidant  flow.  So,  it  is  realistic  to  treat  the  total  heat  capacity  as  a  constant  (assumption  1 1  in  Table  2)  3 

™a(«)cP,a(“)  +'M“)cP,c(n)  =«h,(°)Cp.a(0)  +'M°)<V(0)  (12) 

Substituting  Eqs.  (11)  and  (12)  into  Eq.  (10),  we  obtain  the  required  expression  for  the  local  temperature 


7'(«;vce„)  =  r(0)  + 


x-v« 


Cp.a(°)"c(°)  +C  (0)lWc(0) 


/,i 


(13) 


Note  that  this  expression  for  T(u\ Fcell)  shows  a  simple,  linear  dependence  on  it.  The  term  (y-~  (?eII )/in«  in  Eq.  (13) 
represents  the  total  heat  production  inside  the  fuel  cell  as  a  function  of  the  fuel  utilization.  The  term  <?  (0)ma(0)  + 
cp  c(0)mc(0)  is  the  heat  capacity  of  the  process  flows,  based  on  the  inlet  gas  compositions  and  the  average  temperature  of 
the  fuel  cell  (see  also  the  footnotes  2  and  3). 

Recall  (hat  the  set  (M(jc),Vcen)  has  to  be  calculated  from  BVP(2),  when  only  the  total  fuel  utilization  n,  is  imposed  as  a 
boundary  condition.  Next,  the  total  increase  in  temperature  AT,  defined  by 


Ar=routlcl-riniet  =  7'(W|;vcen)-7-(0) 

follows  from  Eq.  (13),  with  u  —  if.  /mother  possibility,  which  will  be  used  in  the  next  section,  is  to  impose  both  the  total 
fuel  utilization  and  the  total  increase  in  temperature  AT.  Then,  not  only  the  set  (S( jr),V^,ell )  has  to  be  calculated,  but  also  the 
amount  of  oxidant  mc(0)  required  for  cooling.  Rewriting  Eq.  (13)  for  u  =  «,  and  using  the  definition  of  AT,  find 


'K(°) 


(  X  ^cell)^out 

Ar-cp.c(  0) 


Cp.a(O) 

Cp.c(0) 


ma(0)  >  0 


(14) 


with  the  total  current  output  /out  equal  to  lmul.  When  the  oxidant  utilization  is  neglected,  the  cell  voltage  Vcei,  is 
independent  of  the  oxidant  supply.  In  that  case,  the  amount  of  oxidant  required  for  cooling  follows  simply  from  Eq.  (14), 
after  the  cell  voltage  has  been  calculated  for  the  imposed  values  of  AT  and  the  cell  current  /Hut.  When  we  would  take  the 
influence  of  oxidant  utilization  on  the  Ncrnst  loss  into  account,  BVP(2)  and  Eq.  (14)  would  have  to  be  solved 
simultaneously. 


3  For  an  MCFC,  the  change  in  the  total  heat  capacity  if  one  mole  IF  reacts  with  a  half  mole  02  is  about  —7  J/K  (650  °C).  Based  on  this  value,  a  more 
accurate  expression  for  the  total  heat  capacity  is  given  by  ma(0)cp _a(0)  +  mc(0k'p  L.(0)  — 7/in5 /IF.  The  last  term  in  this  relation  represents  the  influence  of 
the  overall  cell  reaction  on  the  total  heat  capacity;  the  influence  of  the  shift  reaction  is  still  neglected.  When  working  with  a  constant  (average)  heat 
capacity,  it  is  recommended  to  use  the  value  that  is  obtained  by  setting  It  =  0.5M]  in  this  relation.  However,  for  MCFCs,  the  difference  with  the  value  given 
by  Eq.  (12)  appears  to  be  relatively  small  (less  :han  1  %  of  the  total  heat  capacity)  and  will  be  neglected  in  order  to  simplify  the  Eqs.  (12)-(14)  as  much  as 
possible.  Nevertheless,  in  both  cases  the  calculus  is  the  same.  Thus  adaption  of  these  equations  is  straightforward. 
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Fig.  3.  Numerical  results  for  the  local  temperatures  inside  an  MCFC  with  co-flow,  obtained  from  a  detailed  numerical  model  [6]  (Tu:  anode  flow.  Tc: 
cathode  flow,  Tc:  electrode-electrolyte  assembly).  For  comparison  also  the  temperature,  7\  as  calculated  from  the  linear  relation  in  Eq.  (13),  is  shown. 


An  analytical  expression  for  the  cell  voltage  will  be  derived  later.  First,  the  accuracy  of  the  simplified  expression  in  Eq. 
(14)  will  be  verified,  using  a  detailed  numerical  cell  model  by  De  Groot  [6].  As  an  example,  we  consider  an  MCFC  with 
standard  gas  supply  (see  Table  3)  at  an  inlet  temperature  of  600  °C.  The  corresponding  values  for  cpa  and  cpc  are  given 
in  Table  4.  In  this  Table  also  the  amount  of  fuel  that  corresponds  with  an  equivalent  input  current  of  1  A  is  given.  From  this 
latter  value  the  necessary  amount  of  fuel  can  easily  be  calculated  if  the  cell  surface,  the  average  current  density  and  the  total 
fuel  utilization  are  known:  for  a  1000  cm2  MCFC,  with  a  current  density  of  157  mA/cm2  (7oul  =  157  A)  and  a  total  fuel 
utilization  of  95%  («,  =0.95),  the  anode:  mass  supply  is  given  by  ma(0)  =  (157/0.95)  X  0.966  X  10“7  =  1.60  X  lO-'7 
kg/s.  Furthermore,  A  T  will  be  set  equal  to  100  K. 

Extensive  numerical  computations  that  take  into  account  variations  in  temperature  over  a  cross  section  of  the  cell  as  well 
as  the  utilization  of  oxidant,  the  heat  release  of  the  shift  reaction  and  changes  in  the  specific  heats  at  constant  pressure  4, 
yield  a  cell  voltage  of  0.75  V.  The  numerically  calculated  value  for  the  oxidant  mass  flow  is  7.02  X  10“4  kg/s,  compared 
with  7.00  X  10“4  kg/s  that  follows  from  Eq.  (14).  For  this  example  the  relative  error  in  the  analytical  expression  for  m/0) 
is  only  0.3%.  In  spite  of  all  the  assumptions  made,  the  oxidant  mass  flow  is  accurately  calculated  from  Eq.  (14),  under  the 
condition  that  the  cell  voltage  is  accurately  known  5 *.  In  the  previous  example,  a  very  high  value  for  the  total  fuel  utilization 
u,  was  imposed  (»,  =0.95).  In  principle,  the  validity  of  Eq.  (14)  is  not  restricted  to  moderate  values  of  but  the 
analytical  expression  for  the  cell  voltage  being  introduced  in  the  next  section  will  not  allow  for  very  high  values  of  ux. 

In  Fig.  3.  the  numerically  calculated  temperature  distribution  is  plotted  as  a  function  of  the  local  fuel  utilization  «(.t).  for 
the  anode  gases  as  well  as  for  the  cathode  gases.  Clearly,  both  temperatures  7a  and  T(  depend  linearly  on  the  fuel  utilization 
in  very  good  approximation.  Solids  and  liquids  are  good  heat  conductors  in  comparison  with  gases.  Consequently,  the 
electrodes  and  the  electrolyte  are  practically  at  the  same  temperature,  say  Te,  at  any  cross  section.  As  indicated  in  the  figure, 
the  local  temperature  of  the  anode  flow  almost  coincides  with  the  hardware  temperature  Te.  It  can  be  shown  that  this  is  due 
to  the  fact  that  only  a  small  amount  of  heat  has  to  be  transferred  to  the  relatively  small  anode  flow.  The  large  cathode  flow 
is  used  for  cooling  and  so  a  serious  temperature  difference  between  the  hardware  and  cathode  gas  will  exist.  Since  the  anode 
gas  temperature  Ta  represents  the  local  temperature  of  the  electrode-electrolyte  assembly  very  well,  it  should  be  used  to 
determine  the  quasi-ohmic  resistance.  However,  the  introduced  mean  gas  temperature  T  is  the  temperature  of  the  gas 
mixture  that  would  be  obtained  by  mixing  of  the  local  anode  and  cathode  flows.  As  illustrated  by  the  figure,  T  is  mainly 
determined  by  the  temperature  of  the  large  cathode  flow. 


4  In  the  numerical  model  [6]  the  specific  heats  cp  a  and  cp  c  are  considered  as  functions  of  both  temperature  and  gas  compositions. 

5  Moreover,  not  all  simplifications  are  essential  in  order  to  obtain  an  expression  for  rnJO).  For  instance,  instead  of  Eq.  (12),  a  more  accurate  average 

value  for  the  total  heat  capacity  can  be  used  (see  Footnote  3).  For  the  MCFC  example,  the  relative  error  in  mc(0)  then  reduces  from  0.3  to  0.1%. 
Furthermore,  in  order  to  eliminate  assumption  7  in  Table  2,  the  last  term  in  Eq.  (14)  can  be  replaced  by  [cpa(0X7^,OP|]ei  —  Tniei)/cp,i:(0XTc  „ulle, — 

Tniet)l"la(0).  Then,  still  one  inlet  temperature  is  assumed  for  both  process  flows  and  A T  is  defined  as  Tt  ouUet  -  7’in,et. 
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3.  Quasi-isothermal  approach  for  an  MCFC  with  co-flow 


3.1.  Problem  definition  for  an  MCFC,  based  on  a  linearized  Nemst  equation  and  a  simple  expression  for  r 


Now  that  we  have  expressed  the  local  temperature  in  terms  of  the  fuel  utilization,  BVP(2)  can  be  reduced  to  a  model 
problem  equivalent  with  the  isothermal  ‘analytical  cell  model’  introduced  in  Ref.  [1],  To  avoid  extended  mathematics,  this 
so-called  quasi-isothermal  approach  will  be  demonstrated  by  a  simple  MCFC  example,  though  the  method  itself  is  quite 
general. 

We  consider  an  MCFC  with  co-flow  and  standard  gas  supply  (see  Table  3).  Of  all  parameters  influencing  the 
quasi-ohmic  resistance  (i.e.  the  polarization  losses),  the  local  temperature  is  the  most  important.  For  reasons  of  simplicity 
we  assume  that  the  quasi-ohmic  resistance  is  a  function  of  the  temperature  only.  Furthermore,  the  temperature  of  the  gas 
flows  increases  from  600  °C  at  the  cell  inlet  («  =  0)  to  700  °C  at  the  outlet  (u  =  tq).  Hence,  the  linear  relation  for  the  local 
temperature,  based  on  the  assumptions  10  and  1 1  in  Table  2,  reads 

r(  m)  =  r^ie,  +  Ar—  =  873  -F  100—  (0  <  if  <  w, )  (15) 

U\  It  | 


Recall  that  the  corresponding  amount  of  coolant  can  easily  be  calculated  once  the  cell  voltage  is  known.  The  purpose  of 
this  section  is  to  obtain  expressions  for  both  the  cell  voltage  and  the  local  current  density. 

The  temperature  dependence  of  the  quasi-ohmic  resistance  r  is  often  described  by  one  or  more  Arrhenius  terms  [10,12]. 
Here,  we  will  focus  on  a  single  Arrhenius  term,  i.e.  we  consider  the  following  simplified  expression  for  the  quasi-ohmic 
resistance 


r(T)  =  r(873)  exp 


1 

873 


(16) 


where  A  is  a  constant.  This  relation  is  convenient  for  a  simple  analysis  and  simulates  actual  quasi-ohmic  resistances  quite 
well,  at  least  for  the  relevant  temperature  range  and  a  suitable  value  for  A.  A  convenient  expression  for  the  reciprocal  value 
of  T,  in  terms  of  n,  is  obtained  by  expanding  Eq.  (15) 


1 

1 

f  100  u 

~  ~  873 

873 1 

1  + - 

I  873  u, 

1 


100  u  \ 

- +0.007  (0  <«<«,) 

873  u  |  J 


(17) 


where  the  relative  error  is  less  than  1%.  Substitution  of  Eq.  (17)  into  Eq.  (16),  yields 


r(  w)  =  r(0)  exp 


(0  <  u  <  tq) 


(18) 


where  A  is  a  constant,  related  to  A.  However,  A  is  fully  determined  if  r(0)  as  well  as  the  value  of  r  at  the  cell  outlet,  i.e. 
r(rq),  are  imposed 


A  = 


( A  >  0;  A  =  0  for  an  isothermal  cell) 


(19) 


In  the  case  of  standard  gas  compositions,  typical  values  are  r(0)  =  1.46  fi  cm2  (600  °C)  and  r(ul)=  0.8  0.  cm2 
(700°C).  Thus,  from  Eq.  (19),  it  follows  that  A  is  typically  0.6. 

Although  the  equilibrium  constant,  K,  for  the  reverse  shift  reaction  is  a  function  of  the  local  temperature,  the  local 
equilibrium  potential  is  reasonably  well  determined  by  using  an  average  value  for  K ,  throughout  the  cell  (i.e.  K  =  0.5  at 
650  °C).  Then  the  expression  for  the  local  equilibrium  potential  in  terms  of  u  and  T  is  the  same  as  for  the  isothermal  model 
(see  Appendix  A  in  Ref.  [1]).  But  now'  the  local  temperature,  T,  depends  on  the  local  fuel  utilization  u  through  the  linear 
relation  in  Eq.  (15) 


RT 

Vju,T]=E(T)-— :/(«;*) 


(20) 


The  explicit  expression  for  the  function  f(u.  K )  is  rather  complicated  and  is  not  presented  here  (see  Ref.  [1]).  The 
reference  potential  E{T)  in  Eq.  (20)  would  be  the  open-circuit  voltage  (OCV)  of  the  cell,  if  T  would  be  the  inlet 
temperature.  The  value  of  E  is  determined  by  thermodynamics  and  appears  to  be  a  linear  function  of  temperature  in  very 
good  approximation  [3],  Therefore,  a  constant  value  for  the  partial  derivative  of  E  with  respect  to  T  will  be  used. 
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In  the  isothermal  model  the  parameter,  a,  was  introduced  to  describe  the  utilization  polarization  (or  Nemst  loss).  It 
represents  a  characteristic  value  for  the  decrease  in  the  local  equilibrium  potential  due  to  fuel  utilization  and  is  defined  as 
the  partial  derivative  of  -  V^[uJ ]  with  respect  to  u;  the  derivative  at  constant  temperature.  Analogous,  the  parameter  a  in 
the  non-isothermal  model  will  represent  a  characteristic  value  for  the  total  derivative  of  -  V^[u.T(u)\  with  respect  to  u. 
Based  on  Eq.  (15).  a  characteristic  value  for  this  total  derivative  derives  from  considerations  made  in  Appendix  A  6 


100/0.5  bE\ 

<5  =  a  H - - a - =  constant 

u !  \  923  'dT  I 


(21) 


Hence,  as  £  is  a  slightly  decreasing  function  of  temperature,  the  parameter  a  in  the  non-isothermal  model  will  be 
somewhat  larger  thanjhe  a  in  the  isothermal  model.  Analogously  to  Eq.  (2),  the  linearized  version  Veq*  of  the  local 
equilibrium  potential  Vcq  in  the  non-isothermal  model  is  given  by 

\/q*  [«,?■(«)]  =Veq*[0,r(0)]  (22) 

Now  that  we  have  expressed  the  local  equilibrium  potential  and  the  quasi-ohmic  resistance  in  terms  of  the  local  fuel 
utilization  u  only,  the  BVP(2)  can  be  solved  for  u  and  Vcell.  To  this  end,  either  the  non-linearized  Nemst  Eq.  (20)  or  the 
linear  approximation  in  Eq.  (22)  can  be  used.  It  is  emphasized  that  numerical  results  will  be  based  on  the  non-linearized 
Nemst  equation,  whereas  analytical  results  will  be  based  on  the  linearized  Nemst  equation.  The  cases  a  =  0  and  a  >  0  will 
be  treated  separately.  The  case  a  =  0  is  hypothetical  and  easy  to  solve.  Its  solution  will  be  used  as  a  reference  solution  to 
show'  clearly  the  effect  of  gas  utilization,  in  practical  cases  where  a  >  0. 


3.2.  Solution  for  a  =  0 


In  the  hypothetical  case  that  the  local  equilibrium  potential  is  independent  of  gas  utilization,  i.e.  a  =  0,  solving  BVP(2) 
is  straightforward.  The  solution  is  then  given  by 


u{  x)  = - ln(  1  -  yA.r), 

A 


Ken  =  Kjo,r(o)]  -<?>/>, 


where  (r)  is  the  average  value  of  the  quasi-ohmic  resistance  and  y  is  a  constant.  Respectively,  they  are  defined  by 
1  rU,  r(0)  —?(«,) 


1  rU  |  i 

(  r)  =  —  /  r(u)  dw  =  — 7— 
u ,  Jo  ln[  r 


y  = 


<r>  l-exp(-A) 
A 


(0)]  -  ln[r(H,)] 


( 0  <  y  <  1 :  y  =  1  for  A  =  0) 


(23) 

(24) 

(25) 


For  an  isothermal  fuel  cell  (A  =  0)  the  constant  y  is  determined  by  the  limiting  value  of  Eq.  (25)  for  A  — >  0. 
An  expression  for  the  local  current  density  easily  follows  from  the  solution  given  by  Eq.  (23) 


du 

'(*)  =  *inT-(*)  = 
djt 


JlB 


1  —  yA.v 


(26) 


with  iout  the  average  current  density,  being  equal  to  imu{.  Since  yA  is  a  constant  smaller  than  unity,  the  local  current  density 
i(.x)  is  monotonously  increasing,  in  agreement  with  the  assumptions  of  a  uniform  equilibrium  potential  and  a  monotonously 
decreasing  quasi-ohmic  resistance.  This  behaviour  of  the  local  current  density  changes  drastically  for  realistic  values  of 
5^0. 


3.3.  Solution  for  a  >  0 

In  general,  BVP(2)  is  rather  complicated  since  the  local  equilibrium  potential  will  be  a  decreasing  function  of  the  fuel 
utilization.  However,  without  a  significant  loss  in  accuracy  the  non-isothermal  model  can  be  simplified  to  an  analytical 
model  that  is  equivalent  to  BVP(l)  for  the  isothermal  model  [1], 


6  In  fact.  K  is  about  0.4  at  600  °C  and  about  0.67  tit  700  °C.  The  inaccuracy  in  Eq.  (2 1 ).  introduced  by  taking  K  =  0.5,  can  be  expressed  by  a  correction 
term  (see  Appendix  A).  For  the  standard  fuel  composition,  it  follows  that  the  a  obtained  from  Eq.  (21)  deviates  considerably  less  than  6  mV  from  the  a 
based  on  the  actual  Kit).  The  effect  of  this  inaccuracy  on  the  final  value  for  the  cell  voltage  will  be  less  than  3  mV  (a  relative  error  of  less  than  0.5  %).  The 
local  equilibrium  potential  will  become  more  sensitive  for  the  precise  value  of  K(T)  if  the  partial  pressure  of  CO  increases.  Hence,  if  the  fuel  is  rich  in 
CO.  the  result  of  Eq.  (21)  may  be  inaccurate. 
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The  quasi-isothermal  problem  formulation  is  based  on  the  mapping  u  ->  u ,  given  by  the  solution  of  the  following  initial 
value  problem 

da  r(  a) 

—  =  (0  <  m  <  m,  )  a(0)=0  (27) 

da  ( r ) 

The  end  value  for  u  is  obtained  by  integration  of  Eq.  (27)  from  a  =  0  to  u  =  a,.  In  this  way  (see  also  Eq.  (24))  it  follows 
that  the  end  values  for  u  and  a  are  ecual.  Hence,  the  mapping  a  -»  u  does  not  change  the  boundary  conditions 

m(0)  =  «(0)  =  0  m(1)=h(1)=«|  (28) 

The  derivative  of  u  with  respect  to  x  is  obtained  from  Eq.  (27)  by  applying  the  chain  rule 

da  da  da  r(  a)  d « 

—  = - =  -7~c —  (29) 

dx  da  dx  <r>  dx 


Eq.  (29)  can  be  used  to  eliminate  the  variable  quasi-ohmic  resistance  and  the  meantime  derivative  du/dx  from  the  ODE 
for  u  in  BVP(2).  In  this  way  the  ODE  for  a  transforms  into  an  ODE  for  a 


—  =  (Kq*  [0,7(0)]  -  5a(«)  -  Vc*,)  (30) 

w'here  Vcel,  is  replaced  by  F.*u  to  indicate  that  the  cell  voltage  is  based  on  the  linearized  Nemst  equation.  By  using  the 
linearized  Nemst  equation,  the  obtained  ODE  can  be  solved  analytically  in  very  good  approximation,  after  which  a  solution 
for  a  is  obtained  by  application  of  the  inverse  mapping  u  -»  u.  For  the  quasi-ohmic  resistance  defined  by  Eq.  (18).  the 
solution  u(ii)  of  Eq.  (27)  and  the  corresponding  inverse  mapping  are  given  by,  respectively 


u(u)  = 


yX 


—  exp  —  A 


u(u)  = - In 

A 


yA- 


Using  the  second  part  of  Eq.  (31),  the  local  equilibrium  potential  in  Eq.  (30)  can  be  expressed  in  terms  of  a 


V  *  [0.7(0)]  -55=  veq*  [0,7(0)]  +  5-^  In 


yA- 


(31) 


(32) 


The  functions  a(x)  and  u(x)  are  equal  at  the  boundaries  (Eq.  (28))  and  will  be  slightly  different  for  intermediate  values  of 
x.  Nevertheless,  consider  also  the  following  linear  approximation  of  Eq.  (32),  obtained  by  setting  ii(x)  =  a(x)  for  all  x 

Kq  [0,7(0)]  -  52  *  veq*  [0,7(0)]  -  5a  (33) 

For  practical  values  of  A  >  0,  this  approximation  for  the  local  equilibrium  potential,  as  a  function  of  the  transformed 
variable  a,  will  be  somewhat  too  sn  ail,  but  nevertheless  reasonably  accurate.  For  instance,  for  A  =  0.6,  5  =  0.2  V  and 
ul  =0.8,  the  maximum  deviation  in  the  approximation  is  easily  calculated  to  be  about  —10  mV,  which  is  small  in 
comparison  with  the  value  of  5. 

In  Appendix  A  it  is  shown  that  substitution  of  Eq.  (33)  into  Eq.  (30)  is  allowed,  if  in  the  latter  equation  the  average  value 
(r)  is  replaced  by  a  somewhat  lower  value  in  order  to  compensate  for  this  artificial  decrease  of  the  local  equilibrium 
potential.  The  replacement  of  <r>  will  be  denoted  by  rm.  In  this  way,  a  problem  formulation  for  a  non-isothermal  fuel  cell 
is  obtained  that  is  much  the  same  as  BVP(l)  for  the  isothermal  model.  Therefore,  this  so-called  quasi-isothermal  problem 
formulation  will  be  referred  to  as  BVP(l'). 


3.3.1.  Quusi-isothermal  problem  formulation 

~  =  —fc  [0,7(0)]  -  5a  -  Vc;„)  (0  <x  <  1) 

d-V  rmlm  V 

a(0)  =0  a(  1)  =  a. 

As  shown  in  Appendix  A,  the  mean  value  rm  for  the  quasi-ohmic  resistance  is  accurately  approximated  by 


1 

f  0  for  A  =  0 

S ) ,  with  S  ~  < 

'  2  2 

i  j  >  0  for  A  >  0 

(34) 

(  5  +  2(7) iin  1 

yA  a 

BVP(  1') 
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The  analytical  expression  for  the  small  correction  term  8  is  fairly  accurate  if  the  artificial  decrease  in  the  local 
equilibrium  potential  (the  difference  between  the  Eqs.  (32)  and  (33))  is  small,  i.e.  for  moderate  values  of  A.  Recall  that 
A  =  0  if  the  quasi-ohmic  resistance  is  independent  of  position,  i.e.  for  the  isothermal  model.  Apparently,  8  is  due  to  the 
non-uniformity  of  r.  In  Appendix  B  it  is  shown  that  a  homogeneous  distribution  of  the  quasi-ohmic  losses  is  favourable 
with  respect  to  the  cell  voltage.  It  is  straightforward  to  verify  that  the  variable  resistance  r  (decreasing  from  cell  inlet  to  cell 
outlet)  will  cause  a  more  homogeneous  distribution  of  the  quasi-ohmic  losses  than  a  uniform  average  resistance  (?)  7 .  In 
Eq.  (34),  the  advantage  of  the  variable  quasi-ohmic  resistance  is  expressed  by  the  term  (1  —  S). 

Note  also  that  the  so-called  quasi-isothermal  approach  is,  in  fact,  a  special  case  of  the  isothermal  problem  formulation  by 
BVP(l)  in  Section  2.1.  That  is,  the  quasi-isothermal  problem  BVP(l')  is  obtained  by  using  the  following  parameter  settings 
in  BVP(l),  with  linearized  Nemst  Eq.  (2) 

r^rm,  a-*  a,  Veq*  (0)  ->  Ve  *  [0,T(0)] 

Once  BVP(1')  is  solved  for  the  set  («(x),Vce*n),  the  actual  u(x )  is  obtained  from  the  inverse  mapping  w(w)  given  by  Eq. 
(31). 


3.3.2.  Solution 

Mathematically  BVP(1')  is  the  same  as  the  isothermal  ‘analytical  cell  model’,  presented  in  Ref.  [1].  Only  the  parameters 
have  a  slightly  different  value,  which  is  indicated  by  the  tilde  ‘  ~  ’.  The  solution  of  the  isothermal  ‘analytical  cell  model’ 
given  in  Ref.  fl],  directly  applies  to  BVP(l')  by  using  the  parameter  settings  given  above. 

For  instance,  a  convenient  and  accurate  expression  for  the  cell  voltage  is  obtained  from  Eq.  (18)  in  [1] 

1  1 

Kell  ~  Kq  [KT(°)]  - -au{  -  rm/uut  -  -Z2rm/oul  (first  order  approximation)  (35) 

where  the  number  Z  is  defined  by 


a 

Z  = - 

2rm*in 

Analogous  to  the  derivation  given  in  the  isothermal  model,  it  can  be  shown  that  the  term  5m,/ 2  in  Eq.  (35)  represents 
the  Nemst  loss.  The  last  two  terms  in  Eq  (35)  are  voltage  losses  due  to  the  quasi-ohmic  resistance  r. 

Next,  consider  the  following  expression  for  ur  obtained  from  Eq.  (10)  in  Ref.  [1]: 


u{x:im)  = 


(0)  -  K*, 


1  —  exp 


—  ax 


The  corresponding  current  distribution  in  the  isothermal  model  follows  from  the  identity  i(x)  =  <jn  du/dx 


,  Kq  [o,r(o)J  -  Ke*n 

=  -  - - : - exP 


a 


(36) 


(37) 


An  expression  for  the  local  current  density  i(x)  in  the  non-isothermal  model  follows  from  Eq.  (36).  by  using  the 
mapping  u(u)  (Eq.  (31)  with  y  =  rm/r( 0))  and  the  identity  i(x )  =  iiBdu/dx  =  i(x)dit/du 


Kx’iin)  ~  '(*An)/ 


r(0) 


-  A 


*OAn) 


(38) 


Eqs.  (35)  and  (38)  for  the  cell  voltage  and  the  local  current  density  are  accurate  as  long  as  the  linearization  of  the  Nemst 
equation  is  accurate,  i.e.  up  to  a  fuel  utilization  of  about  80%.  The  following  examples  will  illustrate  this. 

Fig.  4  shows  the  cell  voltage  of  the  non-isothermal  MCFC,  as  a  function  of  the  total  fuel  utilization  according  to 
several  approximations.  All  the  results  correspond  with  an  average  output  current  density  ioul  equal  to  150  mA/cm2.  Curve 
A  is  obtained  from  the  analytical  expression  in  Eq.  (35).  Numerical  computations  (curve  B)  clearly  show  the  accuracy  of  the 
analytical  results  for  «,  <  0.8.  At  the  current  state-of-the-art,  practical  values  for  w,  usually  lie  between  0.75  and  0.9.  Even 
up  to  an  utilization  of  90%,  the  analytical  results  show  reasonably  good  agreement  with  the  numerical  computations. 
However,  in  practice,  the  cell  voltage  may  deviate  earlier  (i.e.  for  smaller  values  of  ut)  and  faster  from  the  analytical 


7  In  both  cases,  the  local  overpotential  decreases  from  inlet  to  outlet.  The  local  quasi-ohmic  losses 
overpotential  and  the  quasi-ohmic  resistance  from  Eq.  (3). 


rp  are  directly  related  to  the  (square  of  the)  local 
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Fig.  4.  Cell  voltage  of  an  MCFC  as  a  function  of  the  total  fuel  utilization  u,  for  a  constant  average  current  density  of  150  mA/cnr.  Curve  A:  from  Eq. 
(35):  curve  B:  from  numerical  computations,  based  on  BVP(2)  with  non-lincarizcd  Ncrnst  equation;  curves  C  and  D:  from  the  isothermal  model  [1], 
(Standard  gas  supply  (Table  3),  oxidant  utilization  is  neglected,  for  parameter  values  see  Table  5.) 


solution  than  the  numerical  results  suggest.  Because,  at  high  fuel  utilizations,  the  partial  pressure  of  hydrogen  is  very  low 
and  an  additional  concentration  polarization  may  occur  as  a  result  of  slow  mass  transport  rates  in  or  near  the  electrodes. 
Evidently,  the  impact  of  this  effect  also  depends  on  the  average  current  density  /out.  However,  since  we  assumed  r  to  be  an 
explicit  function  of  temperature  only,  the  additional  concentration  polarization  is  not  incorporated  in  our  model  explicitly. 
Anyway,  for  (very)  high  values  of  u{  or  ioul,  the  imposed  values  for  r( 0)  and  Km,),  or  even  the  entire  expression  for  r(u), 
has  to  be  corrected  in  order  to  obtain  a  more  realistic  picture.  But  this  step  is  beyond  the  scope  of  this  article  and  will  be 
omitted. 

In  Fig.  4,  also  the  cell  voltage  according  to  the  isothermal  parameters,  instead  of  the  quasi-isothermal  parameters,  is 
plotted  (curves  C  and  D).  The  parameter  values  for  the  several  curves  are  given  in  Table  5.  The  analytical  value  a  =  0.22  V 
presented  in  this  Table  corresponds  with  a  total  fuel  utilization  of  80%  and  a  negligible  oxidant  utilization.  In  practice, 
oxidant  utilization  will  have  a  small  influence  and  consequently  the  actual  value  of  a  will  be  somewhat  higher  (numerical 
calculations  that  take  into  account  the  utilization  of  oxidant  yield  an  a  of  about  0.24  V).  The  parameter  C  that  appears  in 
Table  5  will  be  explained  in  Section  3.4. 

Fig.  5  shows  the  current  distribution  inside  the  non-isothermal  MCFC  for  a  total  fuel  utilization  of  80%  and  an  average 
current  density  of  150  mA/cm2.  Curve  A!  is  obtained  from  the  analytical  approximation  in  Eq.  (38)  and  shows  reasonably 
good  agreement  with  numerical  computations  (curve  B).  The  deviations  are  mainly  due  to  the  linearization  of  Fq.  (32), 
which  was  the  origin  of  the  artificial  loss  in  the  equilibrium  potential.  Consequently,  the  analytical  approximation  yields  an 
average  current  density  that  is  sometimes  too  small.  Curve  A2  is  obtained  from  a  very  accurate  analytical  expression,  which 
will  be  introduced  next. 


Table  5 

Parameter  values  accompanying  Figs.  4  and  5  (standard  gas  supply,  7X0)  =  600°C) 


Curve 

ocv 

(V) 

a 

(V) 

a 

(V) 

r 

(O  cm2) 

rm 

(It  cm2) 

3£/ar 

(V/K) 

C 

(V) 

A 

lC  [0,7X0)]  =  1.045“ 

0.18 

0.22  b 

1.05d 

-0.02  X  10‘2 

0.003 

B 

Cq[0T(0)]=  1.07 

-0.02  X  10'2 

C 

Vcq*[0;923]  =  1.03 

0.18 

1.08 c 

0.003 

D 

Veq[0;923]=  1.06 

1.08 c 

“  For  h,  =0. 

8.  Determined  from  Eq.  (22)  and  Ihe  condition  Ve  *  [ u |  /2,T(it] 

/2)]=Ve*[n, 

/2;923],  using  Veq[», 

/2;923]  =  0.958  V  as 

a  result  from  Eq. 

(2). 

b  Obtained  from  Eq.  (21)  for  ut  =0.8. 

c  The  quasi-ohmic  resistance  at  650  aC  is  given  by  r  =  ?(«,  /2)  =  1 .08  ( 11  cm"),  from  Eq.  (18). 

d  Value  derived  from  Eq.  (34).  using  A  =  0.6  aid  y  =  0.75  (Eq.  (25)),  corresponding  with  HO)  =  1.46  and  /'(»,)-  0.8  (Cl  cm2). 
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Fig.  5.  Current  distribution  in  a  non-isothermal  MCFC,  according  to  several  approximations.  Curve  A,:  from  Eq.  (38);  curve  A,:  from  Eq.  (40);  curve  B 
shows  the  results  of  numerical  computations,  based  on  BVP(2)  with  non-linearized  Nernst  equation.  The  average  current  density  is  150  mA/cm2  at  a  total 
fuel  utilization  of  80%.  (Standard  gas  supply  (Table  3).  7'(0)  =  600  °C,  AT  =  100  K,  oxidant  utilization  is  neglected,  for  parameter  values  see  Table  5  (A 
and  B)). 


3.4.  Extended  relations  for  the  cell  voltage  and  the  local  current  density 


Inaccuracies  in  the  analytical  expression  for  the  cell  voltage  (Eq.  (35))  occur  as  a  consequence  of  the  so-called  initial  dip 
in  the  non-linearized  Nernst  equation.  This  initial  dip  is  the  deviation  of  the  equilibrium  potential  from  its  linear 
approximation  for  small  fuel  utilizations  and  was  neglected  in  the  analytical  approach.  The  obtained  expression  for  the  local 
current  density  (Eq.  (38))  is  somewhat  less  accurate  than  the  expression  for  the  cell  voltage,  since  it  is  also  based  on  the 
linearization  of  Eq.  (33).  First,  the  inaccuracies  due  to  this  extra  linearization  will  be  diminished.  Next,  a  slightly  modified 
expression  for  the  cell  voltage  will  be  introduced,  that  accounts  for  the  influence  of  the  initial  dip  in  the  Nernst  equation. 

It  is  convenient  (see  Appendix  A)  to  introduce  a  quantity,  /?*',  which  deviates  slightly  from  the  equivalent  input  current 
density  /in 


Z 


Z+8 

The  extended  relation  for  the  local  current  density  is  referred  to  as  ;ext  and  is  given  by 

where  i  is  defined  by  Eq.  (38)  and  T(x;/jnxl)  is  defined  by 


y(xuT)  = 


2  Z 


u1-2u(x,C)]z+. 


+  t  ,n  1  -  yx 


(39) 


(40) 


Note  that  a  times  the  last  factor  in  the  expression  for  T(x)  is  equal  to  the  artificial  loss  in  equilibrium  potential,  i.e.  the 
difference  of  the  Eqs.  (32)  and  (33). 

The  deviation  of  the  extended  relation  i„Kt  from  the  numerical  computations  (see  Fig.  5)  is  mainly  due  to  the  initial  dip  in 
the  non-linearized  Nernst  equation  (analogous  to  Fig.  3  in  Ref.  [1]). 

Finally,  the  influence  of  the  initial  dip  on  the  cell  voltage  is  expressed  by  an  additional  term  C/uv  where  C  is  the 
surface  between  the  initial  dip  and  the  linear  approximation  for  the  equilibrium  potential.  The  extended  expression  of  the 
cell  voltage  as  a  function  of  the  total  fuel  utilization  u,  is  obtained  from  Eq.  (35),  analogously  to  the  derivation  of  Eq.  (22) 
in  Ref.  [l] 


i  . 
+  -z- 

3 


fin  Cut  4”  C 


(0.2  <  h,  <0.8) 


(41) 


Note  that  C  represents  the  same  quantity  as  in  the  isothermal  model,  since  it  is  not  attended  by  a  ‘  ~  ’.  Actually,  the  value 
of  C  is  slightly  influenced  by  local  temperature  gradients.  However,  this  dependence  is  so  small  that  it  has  been  neglected. 
The  correction  term  C/m,  is  especially  significant  for  values  of  u ,  between  0.2  and  0.6. 
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Eq.  (41)  is  a  very  accurate  approximation  for  the  cell  voltage  in  BVP(2)  with  non-linearized  Nemst  equation.  In  fact,  the 
approximation  for  the  local  current  density  in  Eq.  (40)  is  based  on  the  linearized  Nernst  equation.  Nevertheless,  also  this 
latter  expression  was  shown  to  be  an  accurate  approximate  solution  for  the  non-linearized  problem,  for  reasonably  high 
values  of  m,  up  to  about  0.8. 


4.  Conclusions 

The  following  conclusions  have  been  drawn  from  the  present  study: 

1 .  A  noil-isothermal  model  for  fuel  cells  with  co-flow  was  developed,  based  on  the  assumptions  1  to  9  shown  in  the 
Tables  1  and  2.  In  this  article  the  attention  was  focused  on  a  simplified  version  of  this  model,  that  is  also  based  on  the 
additional  assumptions  10  and  11  in  Table  2.  If  these  additional  assumptions  are  also  satisfied,  the  local  temperature  inside 
the  cell  depends  linearly  on  the  local  fuel  utilization. 

2.  The  linear  relation  for  the  temperature  simplifies  considerably  the  non-isothermal  fuel  cell  modeling.  It  is  used  to 
determine  the  analogous  of  the  parameter  a  that  describes  the  Nemst  loss  in  the  isothermal  model.  The  linear  relation  is 
also  suitable  to  express  the  amount  of  oxidant  required  for  cooling  as  a  function  of  the  cell  voltage,  the  cell  current  and  the 
total  increase  in  the  cell  temperature.  The  accuracy  of  this  latter  expression  was  demonstrated  for  an  MCFC,  using  results  of 
extended  numerical  computations.  At  a  total  fuel  utilization  of  95%,  the  deviation  in  the  analytically  calculated  value  for  the 
oxidant  mass  supply  was  shown  to  be  in  the  order  of  0.3%  only. 

3.  A  quasi-isothermal  approach  we.s  introduced.  This  approach  provides  the  possibility  to  generalize  the  results  and 
conclusions  drawn  from  the  isothermal  model  [l],  to  non-isothermal  fuel  cells.  A  quasi-isothermal  parameter  setting  was 
introduced  that  can  be  used  in  the  isothermal  model,  in  order  to  determine  analytical  relations  for  non-isothermal  fuel  cells. 

4.  Analytical  expressions  for  the  local  current  density  and  the  cell  voltage,  derived  with  the  use  of  the  quasi-isothermal 
approach,  were  verified  to  be  very  accurate  up  to  high  fuel  utilizations,  by  comparison  with  numerical  calculations. 

5.  A  close  upper  bound  for  the  cell  voltage  of  a  non-isothermal  fuel  cell  was  derived  and  it  was  shown  that  this  upper 
bound  is  achieved  if  and  only  if  the  quasi-ohmic  losses  are  homogeneously  distributed  over  the  cell  (Appendix  B).  A 
homogeneous  or  almost  homogeneous  distribution  of  the  quasi-ohmic  losses  is  therefore  favourable  with  respect  to  the  cell 
voltage  and  fuel  cell  efficiency. 


5.  List  of  symbols 


c  a  specific  heat  at  constant  pressure  of  the  anode  gas  mixture  J /kg 

cpc  specific  heat  at  constant  pressure  of  the  cathode  gas  mixture  J/kg 

C  surface  between  the  initial  dip  in  V  (u)  and  the  linearized  function  Vt*(u)  V 

E{T)  reference  potential  V 

A  H  enthalpy  change  per  unit  lime  W 

A //thermal  enthalpy  increase,  per  unit  time,  due  to  temperature  rise  W 

A  /^isothermal  reaction  heat  at  constant  pressure  per  unit  time  W 

F  constant  of  Faraday  C/mol 

i  local  current  density  A/m2 

/in  total  current  equivalent  of  the  fuel  supply  A 

iia  total  current  equivalent  of  the  fuel  supply  per  unit  cell  surface  A/m2 

/out  total  current  output  of  a  fuel  cell  A 

iout  total  current  output  of  a  fuel  cell  per  unit  cell  surface;  average  current  density  A/m2 

K  equilibrium  constant  for  concentrations  in  the  reverse  shift  reaction 

anode  mass  flow  kg/s 

mc  cathode  mass  flow  kg/ s 

Q  transfer  of  heat  to  a  system  per  unit  time  W 

r  uniform  quasi-ohmic  resistance  in  the  isothermal  model  fl  nr 

r  non-uniform  quasi-ohmiciesistance  in  the  non-isothermal  model  id  nr 

(?)  average  value  of  the  quasi-ohmic  resistance  r 

O  nr 

rm  a  mean  value  of  the  quasi-ohmic  resistance  r  (see  Eq.  (34))  0  nr 

R  universal  gas  constant  J/mol  K 

T  local  temperature  for  the  gas  flows  K 
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AT  total  temperature  difference  in  a  fuel  cell  (inlet  to  outlet)  K 

u(x )  fuel  utilization  in  the  isothermal  model 

u(x)  fuel  utilization  in  the  non-isothermal  model 

u(u)  invertable  mapping  that  results  in  a  quasi-isothermal  problem  formulation 

u j  total  fueli  utilization 

Veq(i/)  Nemst  potential  in  the  isothermal  model  V 

V  [u,T]  Nemst  potential  in  the  non-isothermal  model  V 

V*(u)  linear  fit  for  V(u)  V 

V'eq’[ii,T(«)]  linear  fit  for  Veq[M,T(u)]  V 

Vcell  cell  voltage  in  the  non-isothermal  model  with  non-linearized  Nemst  equation  V 

Vce*,i  cell  voltage  in  the  non-isothemial  model  with  linearized  Nemst  equation  V 

W  delivered  power  W 

x  scaled,  i.e.  relative,  distance  to  the  cell  inlet 

Z  dimensionless  number 


Greek  letters 


a  characteristic  value  for  the  derivative  of  —  V  (w)  with  respect  to  u  V 

a  analogous  of  a,  for  use  in  the  non-isothermal  model  V 

y  constant  related  to  A  (see  Eq.  (25)) 

A  constant  determined  by  KO)  and  r(ut)  (see  Eq.  (19)) 

actiivation  energy  for  the  conduction  J /mol 

X  constant  of  proportion  in  the  simplified  expression  for  the  reaction  heat  V 


Abbreviations 

MCFC  molten  carbonate  fuel  cell 

BVP  boundary  value  problem 

ODE  ordinary  differential  equation 

Symbols 

~  refers  to  the  non-isothermal  model 

*  attends  potentials  that  are  based  on  a  linearized  Nemst  equation 

ext  extended  relation 
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Appendix  A 

A.l.  The  parameters  5,  x •  rm  an^  l7n' 


A.  1.1.  The  parameter  a,  based  on  the  linear  relation  for  T 

In  the  non-isothermal  model,  the  parameter  — a  represents  a  characteristic  value  for  the  total  derivative  of  the 
equilibrium  potential 


5  =  constant  ~ 


[  du 


[u,T(u)] 


d T  8V„ 


^  r  - 


du  8T 


[h,T(m)]  for  most  u  e  (0,1 ) 


(Al) 


where  the  index  T  indicates  that  the  derivative  holds  for  a  fixed  temperature.  This  derivative  is  equal  to  the  parameter  a  in 
the  isothermal  model  by  definition.  Furthermore,  the  derivative  of  T  with  respect  to  u  easily  follows  from  the  linear  relation 
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for  the  local  temperature  derived  in  Section  2  (Eq.  (13)).  For  instance,  for  the  MCFC  considered  in  Section  3,  based  on  Eq. 
(15),  the  expression  for  a  can  be  rewritten  as 


a  ~  a  — 


100  < 
«,  37 


[“>7X5)] 


(  A2) 


The  partial  derivative  of  the  equilibrium  potential  with  respect  to  7  is  determined  from  Eq.  (20)  (recall  that  7  was 
considered  to  be  independent  of  7,  in  this  equation).  Also  using  the  linear  approximation  for  the  equilibrium  potential  in  the 
isothermal  model,  at  an  average  temperature  of  923  K,  we  successively  obtain 


3V 


3  7 


\uj(u)]  = 

37  v  37 


R  dE  au 

f(u.K)  » - 


27 


37  923 


(A3) 


Analogously  to  the  determination  of  the  parameter  a  in  the  isothermal  MCFC  model  [  I  ],  it  is  convenient  to  determine  a 
in  the  point  of  inflection  of  the  equilibrium  potential.  Since  this  point  must  be  in  the  neighbourhood  of  u  =  0.5  [1],  we  set 
it  =  0.5  in  Eq.  (A3).  The  result  approximates  the  average  value  of  3/,q/37  over  the  interval  [0,  1],  Substitution  of  this 
approximate  value  into  Eq.  (A2)  yields  a  suitable  characteristic  value  for  5  (see  also  Eq.  (21)  Note  that  a  small 
improvement  can  be  achieved  if  /(o,s,k)  is  more  closely  approximated  by  adding  the  difference  between  Veq  and  Veq(0). 

So  far,  the  temperature  dependence  of  K  was  neglected.  For  the  sake  of  completeness  we  will  now  investigate  the 
inaccuracy  due  to  this  simplification.  When  the  K  in  Eq.  (44)  is  considered  as  a  function  of  temperature,  this  equation  can 
still  be  elaborated  as  previously  explained.  The  chain  rule  for  differentiation  then  yields 


-^H“’7(h)] 


37 

37 


au 

923 


d  K  3a  _ 

- u 

d  7  <9  K 


(A4) 


A  typical  value  for  the  derivative  (d7/d7)(3a/37)  is  obtained  from  the  linear  approximation  [a(7  -  0.67)  -  a(7  = 
0.4)]/ 100.  where  a(7 )  follows  from  the  isothermal  model  at  650  °C  by  using  adapted  values  for  K ,  corresponding  with 
700  or  600  °C,  respectively.  For  the  standard  fuel  composition  (Table  3),  the  derivative  is  typically  10~4  X  V/K  and  the 
last  term  in  Eq.  (A4)  contributes  about  6  mV  to  the  value  of  a  (m,  =  0.8). 

Nevertheless,  for  the  standard  fue’  composition,  the  actual  difference  in  the  calculated  values  for  a  appears  to  be 
considerably  less  than  6  mV.  Because,  in  fact,  E  depends  on  7(7)  and  so  the  values  for  37/37  in  Eqs.  (A3)  and  (A4)  are 
different.  Consequently,  for  the  standard  fuel  composition,  these  equations  generate  almost  the  same  value  for  BVcq/hT. 


A.  1.2.  The  parameter  for  the  cell  reaction  H2  +  1  /202  —>  H20 

If  the  enthalpy  of  formation  A/ij^J3  of  1  mole  of  water  is  known,  the  corresponding  value  of  \  can  easily  be  calculated. 
For  instance,  at  600  °C  and  atmospheric  pressure,  A/j^°  is  equal  to  -247  kJ/mol  (based  on  the  JANAF  tables).  The 
electrical  current,  corresponding  with  the  electrochemical  conversion  of  hydrogen  is  given  by  Faraday's  law,  being  2  7 
A/mol  hydrogen.  The  parameter  —  x  was  defined  in  Section  2.2.2,  as  the  ratio  of  the  reaction  heat  and  the  cell  current. 
Hence,  for  the  reaction  of  hydrogen  with  oxygen,  \  's  given  by 


X  =  ~ 


A  h 


H.O 


27 


1.28V 


If  also  the  shift  reaction  occurs,  it  is  convenient  to  use  a  slightly  adapted  value  for  x •  Below,  a  suitable  correction  term 
will  be  introduced  that  depends  on  the  inlet  fuel  composition  and  the  total  fuel  utilization  iq.  The  mole  fractions  of  H2  and 
CO  at  the  anode  inlet,  at  equilibrium,  are  respectively  denoted  by  7,(0)  and  X2(0).  Furthermore,  A refers  to  the 
reaction  heat  of  the  shift  reaction  (  —  36  kJ/mol  at  600  °C  and  atmospheric  pressure).  In  order  to  calculate  the  amount  of 
CO  converted  by  the  shift  reaction  we  consider  two  alternative  subprocesses.  In  the  first  subprocess,  the  inlet  fuel 
composition  is  heated  to  the  outlet  temperature  of  the  fuel  cell  and  the  gas  composition  changes  only  as  a  result  of  variations 
in  7(7).  The  amount  of  moles  of  CO  converted  in  this  first  step  is  given  by  72(0;7inlel)  —  X2(0;  70U(lct)  times  the  total  (and 
constant)  number  of  moles  in  the  fuel  flow.  The  constants  7inlcl  and  70UlleI  represent  the  equilibrium  constant  for  the 
reverse  shift  reaction  at  the  inlet  and  outlet  temperature,  respectively.  In  the  second  step,  the  fuel  is  converted  by  an 
isothermal  fuel  cell.  The  amount  of  CO  converted  in  this  latter  step  is  proportional  to  X2(0;7olltlet)  X  v(u] :  7<jllllct)  (the 
explicit  expression  for  the  reaction  coordinate  u(m,)  is  given  in  Ref.  [1]).  The  modified  value  of  x  represents  the  ratio  of 
the  total  reaction  heat  and  the  cell  current 


A/tj^  I  X2(0;7mlel)  -X2(0;7outlel)[l  -r(u,;70U„e,)]  )  A/t^' 
27  \  [X,(0)+X2(0)]M|  j  IF 
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A.  1.3.  The  parameter  rm 

In  fact,  the  mean  value  rm  for  the  quasi-ohmic  resistance  is  the  limiting  value  of  the  following  combined  iteration 
process 


d  u 
d  x 


1 


i  =  — (Kq*[°;n°)]  -Vce*„),  «,-(<))  =0,  i  =  3,4, 5,... 

-  an,-  VJ„ 

d w,,  r  =  3,4,5, .. . 


/■».  ve;[o;no)]  -g«, 

r'  +  l”  «>  '<>  V-  [0;T(0)]  yc* 


with  the  following  initial  values,  defined  in  Section  3 


u2(x)  =  u(x)  (0<,v<l)  r3  =  (r) 

Recall  that  the  constant  «,  is  used  to  denote  the  total  fuel  utilization.  The  derivation  of  the  ODE  for  i  =  3  has  been 
explained  in  Section  3,  omitting  the  index  ‘3’.  It  was  shown  how  the  variable  quasi-ohmic  resistance  in  the  ODE  for  u  is 
replaced  by  its  average  value  r3  =  <r>.  Attention  was  also  paid  to  the  linearization  of  the  ODE  for  i  =  3  and  the 
corresponding  artificial  loss  in  equilibrium  potential.  All  next  iterations  are  based  on  the  same  principle.  However,  now  an 
exact  correction  for  the  artificial  loss  in  equilibrium  potential  is  made,  by  a  local  adaption  of  the  average  quasi-ohmic 
resistance,  instead  of  a  global  adaption.  This  yields  a  ‘quasi-ohmic  resistance’  that  depends  on  position  again,  which  give 
rise  to  the  next  iteration  cycle:  the  calculadon  of  the  average  value  r4  and  the  function  «4. 

Theoretically  this  process  can  be  continued,  and  the  quantities  u3(x).  r4,  «4(.v),  r5,...  can  successively  be  obtained 
from  the  iteration  process.  For  the  same  reason  as  given  in  Section  3  for  the  first  step,  the  end  condition  «,(1)  =  «,  will  hold 
for  all  i.  Consequently,  all  r,  must  lie  within  limits  that  are  independent  of  i.  Because  otherwise  the  end  condition  cannot 
be  satisfied  for  some  z,  since  the  derivative  dzq/d.t  in  the  corresponding  ODE  would  be  too  large  (or  too  small). 
Furthermore,  it  can  be  shown  that  ;<,(x)  >  for  all  i,  and  that  r;  is  monotone  decreasing  with  respect  to  i  8.  So,  r: 

converges  to  a  limit  rm  as  i  tends  to  infinity.  This  implies  that  w,  converges  to  a  limiting  function,  which  satisfy  BVP(l')  in 
Section  3.  This  limiting  function  will  be  denoted  by  ux.  In  Section  3.3,  an  accurate  approximation  for  rm  was  used,  which 
will  be  derived  below.  The  expression  for  the  cell  voltage  derived  in  Section  3.3  is  accurate.  But  the  difference  between 
ti3(.v),  referred  to  as  u(x),  and  ux(x)  was  neglected.  Consequently,  the  expression  for  the  local  current  density,  derived  in 
Section  3.3,  is  somewhat  less  accurate,  since  it  is,  in  fact,  based  on  the  function  ux  instead  of  the  required  u. 

Assuming  that  the  maximum  difference  between  u3(x)  and  ux(x )  is  small,  an  accurate  analytical  approximation  for  rm 
can  be  derived  as  follows.  Dividing  the  ODEs  for  ujx)  and  n3(x)  we  obtain  (also  using  the  initial  settings  for  «,  and  r3) 

d«x  =  < r>  ^[0:^(0)]  -  a “x- Kelt 

d«3  /m  Ve*[0;T(0)]  -55 


By  setting  ujx)  =  «3(.r)  on  the  right-hand  side  of  this  equation  (recalling  that  uj.x)  >  w3(v))  and  integrating  the  result 
from  m3  =  0  to  h3  =  we  obtain 


r 


m 


Kq*[o;r(o)]  -a«3-v«*u 

— - — — d  U 

Ve*[0;r(0)]  -««(«,)  -Vce*„ 


( A5) 


where  also  the  end  condition  uS  1)  =  u{  is  used.  Next,  from  a  suitable  standard  expansion  of  the  right-hand  side  of  this 
equation,  it  can  be  shown  that  the  following  upper  bound  is  fairly  close  to  the  exact  value  of  rm 


1 


&  rU\  _ 

—  I  [  u  —  u(  u)  jd« 
M)  •'() 


(A6) 


with 


a  2  a 

^To;r(0)]  -  yc;u  ~  (5  +  2<r)«,..)«I 

The  approximation  for  a  is  based  on  the  so-called  double  linear  relation  or  zero  order  approximation  for  the  cell  voltage, 
obtained  by  omitting  the  last  term  in  Eq.  (35)  (see  Ref.  [1]:  without  a  further  loss  in  accuracy  rm  is  replaced  by  <r)  in  the 
approximation  for  a). 


8  Under  the  condition  that  r(u)  is  a  monotone  decreasing  function  of  u.  In  the  general  case,  the  inequality 
proof  of  convergence  is  somewhat  more  extended. 


is  not  necessarily  valid  and  the 
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In  specific  cases,  the  expression  for  rm  given  in  Eq.  (A6)  (or  a  more  complete  expansion  of  Eq.  (A5))  can  be  elaborated 
by  using  the  inverse  mapping  u(u).  The  expression  for  rm  that  was  given  in  Section  3.3  (Eq.  (34)).  is  obtained  by  using  the 
mapping  in  Eq.  (31). 


A.  1.4.  The  parameter 

The  extended  expression  for  the  local  current  density  given  in  Section  3.4  is  obtained  from  a  linearization  around  the 
limiting  solution  ux.  This  linearization  makes  it  possible  to  bypass  the  iterations  for  /  >  4  in  a  more  accurate  way  than  by 
neglecting  the  difference  between  w,  and  ux.  Here,  the  main  point  of  this  exercise  will  be  highlighted,  i.e.  the  introduction 
of  the  parameter 

The  limiting  solution  ux  is  based  on  the  parameter  value  rm,  whereas  the  function  «,  is  based  on  the  parameter  value 
(r).  However,  both  functions  correspond  with  the  same  cell  voltage,  since  this  quantity  remains  unchanged  by  the  iteration 
process.  Hence,  both  functions  must  also  correspond  with  the  same  Total  resistance’  that  determines  the  cell  voltage.  In  the 
isothermal  model  (see  Eq.  (25)  in  Ref.  [1])  it  was  shown  that  this  total  resistance  consist  of  the  load  resistance,  the  uniform 
quasi-ohmic  resistance  and  the  ‘utilization  resistance’,  being  equal  to  a/2im.  The  utilization  resistance  depends  on  the 
equivalent  input  current  density  ;in.  V/hen  switching  from  a  quasi-ohmic  resistance  (r)  to  a  quasi-ohmic  resistance  rm  (i.e. 
bypassing  the  iterations  for  i  >  4)  the  equivalent  input  current  density  has  to  be  adapted  in  such  a  way,  that  the  total 
resistance  remains  constant.  Of  course:  this  is  a  mathematical  trick  and  does  not  affect  the  real  value  of  /jn.  So,  the  parameter 
easily  follows  from  the  required  change  in  the  utilization  resistance,  determined  by 


a 


+  (r)  = 


2  if: 


+  r„ 


(A7) 


It  can  be  made  clear  that  this  equaiion  for  /inext  has  to  be  be  used  in  combination  with  Eq.  (A6)  for  rm.  Recall  from  Eq. 
(34),  that  the  parameter  rm  is  equal  :o  the  average  resistance  (r)  times  the  factor  (1  -  5).  When  rewriting  Eq.  (48),  one 
may  approximate  5/(1  —  5)  by  8  in  order  to  obtain  Eq.  (39). 

Eq.  (A6)  gives  an  approximation  for  rm,  hence  for  8.  The  obtained  expression  for  i“*  is  recommended  if  this 
approximation  for  8  is  used.  However,  if  8  is  more  accurately  determined,  it  is  recommendable  to  use  the  identity 
t'fnxt  =  fin/v  1  -  8  instead  of  Eq.  (39).  It  can  be  shown  that  this  latter  identity  yields  the  most  accurate  result,  as  long  as  5  is 
very  accurately  known. 


Appendix  B 

B.l.  A  close  upper  bound  for  the  cel!  voltage  in  BVP(2) 

Denoting  the  average  value  of  the  square  root  of  the  quasi-ohmic  resistance  by  (\f?)  (analogous  to  Eq.  (24)).  it  will  be 
shown  that 


Kell  ^  /  Kq[«,rOl)]d«  -  <v^>2r 

«1  Jo 


(A8) 


with  equality  if  and  only  if  r(x)r(x)  =  constant  for  all  x.  The  upper  bound  in  Eq.  (A8)  is  the  generalization  of  the  upper 
bound  for  the  cell  voltage  that  was  derived  from  the  isothermal  model,  hence  for  the  special  case  in  which  r(x)  =  r  for  all 
A'  [1]. 

In  order  to  show  the  validity  of  Eq.  (A8)  we  multiply  both  sides  of  the  ODE  in  BVP(2)  by  du/dx  and  integrate  with 
respect  to  x  from  zero  to  unity.  Rewriting  the  result,  using  Eq.  (4).  we  obtain 


Kell  ~  ~  f  'v  [3  J(B)]du~—  fir(x)P(x)d: 

«1  A)  'out  -A) 


(A9) 


The  upper  bound  in  Eq.  (A8)  is  obtained  by  estimating  the  last  integral  in  Eq.  (A9). 

This  estimation  is  based  on  the  Lemma  of  Cauchy-Schwarz  [18,19],  which  states  that  the  following  inequality  holds  for 
two  continuous  functions  f(u)  and  g(t7) 


v  2 

f  '/(  u)g(u)du  |  <  f  \f2(u)du(  'g2(u)du 

Jet  I  J() 


(  A10) 
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Applying  the  Lemma  of  Cauchy-Schwarz  on  the  following  integral 


by  setting  /(«)  =  y 'r(fi)  y 'i(u)  and  g(u)=  1  /\!i(u)  ,  we  obtain 
(if?)2  <—(f  r(u)i(u)du  f  — - d«j 

1/0  J()  l(u)  J 

Using  the  identity  i(x)  —  ijndu/dx  (Ec|.  (4)),  this  inequality  can  be  rewritten  as 


(All) 


(A12) 


Substitution  of  this  inequality  into  Eq.  (A9)  yields  Eq.  (A8). 

The  Lemma  of  Cauchy-Schwarz  also  states  [19],  that  the  equality  in  Eq.  (A  10)  holds  if  and  only  if  /=  constant  X  g. 

Hence,  the  equality  in  Eq.  (All)  holds  if  and  only  if  =  constant/ v'f.  So,  the  cell  voltage  is  maximum  if  the 
quasi-ohmic  losses  are  homogeneously  distributed  over  the  cell 

r(x)r(  x)  =  constant  for  all  x  <=>  the  upper  bound  for  Vcell  inEq.  (A8)  holds 


For  the  MCFC  considered  in  Section  3.  and  a  total  fuel  utilization  of  80%,  the  deviation  of  the  cell  voltage  from  the 
obtained  upper  bound  is  10  mV  (from  Eq.  (41)). 
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